In the following notebook, the sound field of an ideal plane wave impinging on a spherical microphone array is generated and visualized. Then, the impulse response of two different looking direction is extraced using plane wave decomposition.
import sys
sys.path.insert(0, '../')
from sound_field_analysis import gen, plot, utils, process
from plotly.offline import init_notebook_mode
init_notebook_mode()
First the limit for the spherical domain is set. Lower order correspond to faster calculation time at a loss of fidelity in the radial spatial domain.
Then, the configuration of the microphone array that is to be simulated. It simply holds three elements:
Furthermore, the wavetype is set to 'plane' and the direction of the wave is specified to come from "straight ahead". NFFT denotes the number of the FFT bins and fs the sampling frequency used for the simulation.
order = 5
array_configuration = [0.5, 'rigid', 'cardioid']
wavetype = 'plane'
azimuth = utils.deg2rad(0)
colatitude = utils.deg2rad(90)
NFFT = 128
fs = 48000
gen.ideal_wave() returns the spherical fourier coefficiens of a ideal plane wave as captured by an ideal array of the provided configuration. To further work with these, it is often helpful to also generate the corresponding radial filters along side.
field_coefficients = gen.ideal_wave(order, fs, azimuth, colatitude, array_configuration, NFFT=NFFT)
radial_filter = gen.radial_filter_fullspec(order, NFFT, fs, array_configuration)
To visualize the field, the spherical fourier coefficients and the generated radial filters are passed to plot.makeMTX() along with a kr bin. The function returns the sound pressure at a resolution of 1 degree, which then can be visualized using plot.plot3D(). Several styles ('shape', 'sphere', 'flat') exist.
kr_IDX = 64
vizMTX = plot.makeMTX(field_coefficients, radial_filter, kr_IDX)
plot.plot3D(vizMTX, style='shape')
In order to extract the impulse response, a plane wave decomposition (PWD) is performed using he process.plane_wave_decomp() function, whos output in the frequency domain can then simply be transformed to the time domain using process.iFFT(). We can calculate several looking directions at once by supplying a vector of azimuth / colatitude pairs.
Since this is an approximated ideal plane wave, the impulse response will be almost dirac impulse-like in the direction of the plane wave (azimuth = 0°, colatitude = 90°) and almost zero at a 90° angle (like azimuth = 90°, colatitude = 90°).
looking_direction = [[0, utils.deg2rad(90)],
[utils.deg2rad(90), utils.deg2rad(90)]]
Y = process.plane_wave_decomp(order, looking_direction, field_coefficients, radial_filter)
impulseResponses = process.iFFT(Y)
spectrum, f = process.FFT([impulseResponses, fs])
plot.plot2D(impulseResponses, type='time', fs=fs)
plot.plot2D(utils.db(spectrum), type='logFFT', fs=fs)